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THE  CUBIC  SMOOTHING  SPLINE  AS  A  DIGITAL  FILTER 


Kenneth  Peters  and  Edward  R.  Cook 


ABSTRACT 

Reinsch  (1967)  introduced  the  cubic  smoothing  spline  as  the  solu¬ 
tion  to  a  certain  variational  problem  in  the  theory  of  data  smoothing. 

This  smoothing  operator  has  good  transient  response,  a  minimal  curva¬ 
ture  property,  and  compensation  for  end  effects  so  that  all  the  data 
is  used.  In  addition,  the  data  can  have  arbitrary  spacing  and  weighting. 

In  this  paper  we  represent  Reinsch1 s  equations  in  terms  of  convolu¬ 
tions  for  the  special  case  of  equally-spaced  and  equally-weighted  data. 

We  then  use  a  Fourier  transform  to  obtain  a  frequency  domain  represen¬ 
tation,  solve  the  characteristic  polynomial  of  the  system,  and  apply  an 
inverse  transform  to  the  frequency  response  function  to  obtain  the 
impulse  response  function.  The  Lagrangian  multiplier  of  the  variational 
problem  parameterizes  the  family  of  digital  filters  derived  in  this  man¬ 
ner.  This  representation  which  is  exact  for  finite  data  sets  extends  the 
range  of  application  of  the  smoothing  spline  and,  incidentally,  simplifies 
the  algebra  and  reduces  the  amount  of  computer  time  and  storage  required, 
by  allowing  one  to  select  the  degree  of  smoothing  on  the  basis  of  frequency 
domain  considerations.  One  can  do  time  domain  filtering  without  end 
effects  in  a  computationally  efficient  manner. 

Smoothings  of  four  tree-ring  width  sequences  from  forest  interior 
sites  are  presented  to  illustrate  the  particular  suitability  of  these 
filters  for  data  exhibiting  both  non-stationary  and  episodic  behavior. 
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THE  CUBIC  SMOOTHING  SPLINE  AS  A  DIGITAL  FILTER 


I  INTRODUCTION 


Reinsch  [7]  discusses  the  problem  of  smoothing  noisy  data  when  the 
underlying  function  is  not  known.  We  begin  with  a  brief  statement  of 
his  main  results  in  his  own  notation  and  terminology.  He  suggests  that 
the  fitting  function,  g(x),  minimize  the  total  square  curvature, 


2 

dx 


(la) 


xonder  the  constraint  that  a  weighted  sum  of  squares  of  the  residuals 
be  bounded  above , 


n 


"g(Xi)-yil  2  <  S; 


6Yi 


(lb) 


where  the  y^  represent  the  noisy  data,  the  6y^  are  a  set  of  weights 
and  S  is  a  scaling  parameter.  Equivalently,  introducing  a  Lagrangian 
multiplier,  p,  and  an  auxiliary  variable,  z,  the  functional 


is  to  be  minimized. 

The  solution  is  a  cubic  spline,  i.e.,  a  concatenation  of  cubic 
segments : 

g^(x)  =  ai  +  (x-x^)  +  c^  (x-x^) 2  +  d^  (x-x^)  3,  x^  ^  x  <  x^+1,  (2) 

such  that  g(x) ,  dg(x)/dx,  and  d2g(x)/dx-  are  continuous  at  x^  for 
i=l,  n-1. 
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If  the  Lagrangian  multiplier,  p,  is  known,  the  values  can  be 
computed  from  ^ 

(QD  Q  +  pT )  c  =  pQy,  (3) 

where  the  circumflex  denotes  transposition,  and  the  a^  from 


Then 


a  =  y  - 


p-1D2Qc. 


(4) 


di  =  (ci+1  -  ci)  /  3hi,  i  =  0,  n  -  1  (5) 

where,  for  the  sake  of  uniformity  of  notation,  Reinsch  chooses  to  set 
Cq  and  c^  equal  to  0 ;  and 


bi  =  E(ai+1  ~  a±) /hi3  -  cihi  -  dihi  ,  i  =  0,  n  -  1. 
The  notation  used  is  Reinsch' s: 


hi  =  * 


i+1  "  xi; 


c  (c]_,  •  •  •  'cn-l^  ' 


(6) 

(7a) 
(7b) 
(7c) 
( 7d) 
(7e) 
(7f) 


y  =  (y0#  yi . yn.r  yn)  * 

a  ^ a0  '  al'  •  •  •  '^-1'  an^  ' 

D  =  diag  (6yQ, . . . ,6yn) ; 

T  =  a  positive-definite  tridiagonal  matrix  of  order  n-1, 

where  t^  =  2  (hi_1  +  hi)/3,  and  ti>i+1  =  h±/3  =  ti+1/i; 

Q  =  a  tridiagonal  matrix  with  n+1  rows  and  n-1  columns  (7g) 

where  qi_1^i  =  l/hi_1,  qi;j_  =  -  (l/hj^  +  l/hi)  ,  and 

gi+l,i  "  1/bi* 

For  Reinsch1 s  intended  application  to  data  smoothing  it  is  natural 
to  specify  the  weights  5y^,  usually  as  the  square  root  of  the  error 
variance  at  x=x^  if  it  is  known  or  can  be  estimated,  and  S,  as  the  value 


of  a  chi-square  variable  with  n  degrees  of  freedom.  He  then  proceeds  to 
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and  may  be  of  interest  to  the  reader.  Reinsch  t  8]  recasts  the  problem  of 
reference  £71  in  much  more  general  terms  and  shows  that  similar  results  hold. 
Demmler  and  Reinsch  £41  achieve  by  advanced  methods  a  very  general  result  of 
which  formula  (16)  of  this  paper  is  a  special  case.  Wahba  £9]  discusses  the 
problem  of  optimizing  the  degree  of  smoothing  given  certain  assumptions 
regarding  the  random  component  of  the  data  and  gives  explicit  formulas 
which  can  be  used  in  conjunction  with  the  heuristic  considerations  which 
motivated  this  paper.  Wahba  and  Wold  £10]  and  Craven  and  Wahba  [3]  apply  a 
method  of  cross-validation  to  select  the  smoothing  parameter  and  in  £11] 
Wahba  and  Wold  discuss  the  use  of  periodic  splines  in  spectral  density 
estimation  which  is  related  to  a  suggestion  for  natural  splines  in  this 
paper  in  Section  V.  Wold  £12]  is  concerned  with  the  suitability  of  splines 
of  various  kinds  from  a  data  analytic  point  of  view.  In  this  paper  we  are 
concerned  with  a  discrete  problem  —  the  relationship  between  the  vectors  y 
and  a.  In  the  limit  of  large  p  there  are  no  digital  filtering  effects. 
Horowitz  C6]  considers  the  continuous  problem  of  the  filtering  effects  of 
common  interpolating  splines.  His  analysis  may  be  of  interest  if  the  dis¬ 
crete  data  arise  from  accumulation  or  from  sampling  continuous  records. 

The  references  given  here  should  provide  an  adequate  introduction  to 
smoothing  splines. 

The  rest  of  the  paper  consists  of  four  sections  and  two  appendices. 
Section  II  contains  the  derivation  of  the  convolutional  representation 
and  of  the  frequency  response  function.  Section  III  contains  the  solution 
of  the  characteristic  equation  of  the  system.  Section  IV  contains  the  deri¬ 
vation  of  the  impulse  response  function.  Section  V  includes  a  discussion 
of  how  to  interpret  and  apply  the  representation  derived  in  sections  II, 

III  and  IV,  and  an  example  which  illustrates  the  time  domain  behavior  of 
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the  cubic  smoothing  spline.  Appendix  A  is  a  discussion  of  this  behavior  in 


the  limit  of  small  curvature,  and  Appendix  B  contains  an  analysis  of  the 
coupled  ordinary  differential  equations  of  which  the  cubic  smoothing  spline 
is  an  interpolated  finite  difference  solution. 

The  reader  interested  only  in  practical  applications  should  read 
sections  I  and  V  and  consider  the  main  practical  results  which  are  in 
equations  (16)  and  (33)  and  figures  1  and  5.  The  critical  reader  will 
want  to  include  the  beginning  of  section  II.  The  appendices  and  the 
parts  of  the  text  relating  to  the  correction  terms  (indicated  by  primes) 
and  the  small  curvature  limit,  can  be  ignored,  but  clues  to  understanding 
the  Fourier  representation  are  contained  therein.  Finally,  those  filters 
generated  by  the  real  negative  roots  of  the  characteristic  polynomial 
are  of  little  practical  interest  since  they  result  in  very  little 
smoothing. 

II  THE  CONVOLUTION  REPRESENTATION  AND  THE  FREQUENCY  RESPONSE  FUNCTIONS 

Consider  the  special  case  of  equations  (3)  and  (4)  for  which 
h^  =  1.0  =  Sy^  for  all  i.  In  this  case  the  matrices  Q  and  T  and  R  = 

QQ  +  pT  will  be  Toeplitz,  i.e.  ,  all  the  terms  on  the  same  diagonal  will 
be  equal,  and  symmetric.  The  matrix- vector  product  Rc  will  be  equivalent 
to  convolution  with  a  symmetric  (phase-free)  filter  except  for  the  first 
and  the  last  two  points.  In  fact  these  restrictions  are  necessary 
because  convolution  is  defined  as  a  time-invariant  linear  process,  and 
it  is  impractical  to  do  Fourier  transforms  on  unequally  spaced  data. 

The  particular  numerical  values  of  the  spacing,  in  this  case  unity, 

-3 

serves  to  scale  or  normalize  the  parameter  p  which  has  the  units  of  x 
With  these  restrictions  then  the  desired  representation  can  be  achieved 
as  follows. 
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First,  extend  a  typical  row  of  the  matrix  R  by  appending  an  infinite 
number  of  zeroes  on  each  end,  and  do  the  same  for  Q  and  Q.  This  exten¬ 
sion  is  necessary  because  the  limits  on  the  Fourier  integral  are  infinite. 
The  particular  form  of  the  extension  (zeroes)  is  dictated  by  the  band 
form  of  the  corresponding  matrices.  Thus  define 

r  =  {...,0,1, (-4+p/3)  ,  (6+4p/3)  ,  (-4+p/3) ,1,0, . . . }  (8) 

and 

q  =  {... ,0,1, -2, 1,0,.. .}  (9) 

in  accordance  with  (7f)  and  (7g) .  The  subscript  p  makes  the  para¬ 
meterization  explicit.  The  vectors  c,  y  and  a  can  be  extended  in  a 
similar  manner,  i.e., 

c  =  {.. . ,0,cQ,c1,... ,cn_1,cn,0, .. .},  c0=cn=0,  (7b') 

y  =  {...,0,y0,y1,...,yn_1,yn,0,...},  (7c1) 

and 

a=  { . . . ,0,aQ,a1, . . . ,an_1,an,0, . . . };  (/a') 

but  this  particular  extension  for  c,  y  and  a  is  not  dictated  by 

Reinsch's  equations  (except  for  the  end  conditions  Cg=cn=0)  which 

entail  no  assumptions  regarding  the  behavior  of  y  outside  the  data 

interval.  We  use  it  only  for  definiteness  and  convenience.  More 

general  extensions  will  be  discussed  in  section  V.  The  main  point 

here  is  to  represent  Reinsch's  computations  as  exactly  as  possible. 

Next  notice  that  while  (3)  equates  two  vectors  of  length  n-1, 

r  *  c,  where  *  denotes  convolution,  and  q  *  y  will  in  general  have  n  +  3 
J? 

non-zero  elements.  To  represent  (3)  accurately  the  extra  terms  gener¬ 
ated  by  the  convolutions  must  be  subtracted  from  both  sides  before  they 
are  set  equal.  Thus,  (3)  becomes 

rp*c  -  =  p  (q  *  y-y')  ( 3 ’ ) 
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with 


c  '  =  {  .  .  .  ,0,0^  ,  i  (-4+p/3)  c^+c^  ,0  , .  . .  ,0,  Ccn_2+  (-4+p/3)  cn_12  , 
cn_]_,0, .  . .  }  (10) 

and 


y'  =  { • • • /0,y0, (-2y0+y1) ,0, . . . ,0, (-2yn_1+yn) ,yn,0, . . . }.  (11) 

Note  that  the  prime  will  not  be  used  to  denote  differentiation  in  this 
paper. 

In  the  case  of  equation  (4)  the  extra  terms  are  retained.  It 
becomes 


y  -  —  q 
P 


c,  p  >  0, 


(4') 


Finally,  to  achieve  symmetry  in  the  transforms,  and  for  representa¬ 
tional  purposes,  the  discrete  vectors  will  be  represented  as  weighted 
sums  of  delta  functions.  Letting  a  denote  any  of  them  define 


a(t)  7  a.  6(t-i).  (12) 

*• — l 

i=-°° 

The  explicit  dependence  on  t  will  identify  continuous  time  domain 
functions. 


Fourier  transforming  (3')  and  (4')  modified  according  to  (12) 


gives 

r  (f)c(f)  -  Cp  ( f )  =  p  Cq(f)y(f)  -  y'(f)] 

(13) 

and 

a (f )  =  y(f)  -  |q(f)c(f)  ,  p  >  0. 

ir 

(14) 

Solving 

(13) 

for  c(f),  substituting  in  (14)  and  rearranging  gives 

a  (f) 

=  u  (f)y(f)  +  u  '  (f)  Cy'  (f)  -  ic  •  (f)]  ,  p  >  0, 
p  P  P  P 

(15) 

where 

u  (f) 

=  (p/6)  cos2TTf+ (p/3) 

(16) 

P 

cos^2'n'f+  (  (p/6)  -2)  cos2nf+  (  (p/3)  +1) 

and 

n  ,  _  2  (cos2TTf-l) 

(17) 

cos^2iTf+  (  (p/6)  -2)  cos27rf+  (  (p/3)  +1 

which  are  the 

frequency  response  functions. 

The 

functions  u  (f)  and  u  '(f)  are  shown  in  figures  1  'and  2  over 
P  P 

the 

useful 

range 

of  p. 
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Ill  THE  ROOTS  OF  THE  CHARACTERISTIC  POLYNOMIAL 


The  characteristic  polynomial  of  the  system  is  the  (common)  denomi¬ 
nator  of  u  (f)  or  u  '  (f)  expressed  in  terms  of  z  =  exp(27rif)  .  The  roots 
I?  P 

of  this  function  will  be  needed  in  section  IV.  Setting  the  denominator 
of  Up(f)  equal  to  zero  and  solving  for  cos2irf  gives 


and 


Zp  =  1-  (p/12)  +  i(p(72-p))1/2/12 


z"  =  1- (p/12)  -  i (p (72-p) ) 1/2/12 


(18) 


(19) 


The  principle  root  is  used  here  and  always.  These  roots  lie  on  the 

circle  of  radius  3  centered  on  the  real  axis  at  -2  for  0<p$72. 

For  p£72  Im  z  =0.  Also,  to  avoid  ambiguity,  lim  z  as  p-*30  is  defined 
P  P 

to  be  -2. 

If  z=exp  (27rif )  is  a  root,  so  is  1/z  since  cos  2nf  =  (exp(2irif)  + 
exp  (-27rif )  ) /2 .  Using  this  relationship  and  its  converse,  the  four  roots 
are,  for  0<p<72, 


and 


++  +  •  +2»  1/2 

z  =  z  +  i  ( 1-  z  ) 

P  P  P 


H —  +  .  , +2 .  1/ 2 

z  =  z  -  i(l-z  ) 

p  p  p 


-+  -  -2  1/2 
Zp  =  zp  +  i(l-zp  ) 


-2  1/2 

z  =  z  -  i(lz  ) 

P  P  P 


(20) 

(21) 

(22) 

(23) 


For  p^72  there  are  four  distinct  negative  real  roots.  Two  of 

these,  z++  and  z  +,  remain  within  the  unit  circle, 

P  P 


and 


z  + 
P 


(24) 

(25) 
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while  the  reciprocals,  z^+  and  z ^  ,  respectively,  remain  outside.  The 

lim  z  as  p-*°°  is  defined  to  be  0,  while  lim  z  as  p-*»  is  -2  +  /3. 

P  P 

The  root  z^"*"  is  plotted  in  figure  3.  The  other  roots  are  simply 
related  to  z^+.  Thus  z^  is  the  reciprocal,  z^~  is  the  complex  con¬ 
jugate,  and  z  is  the  reciprocal  of  the  complex  conjugate. 

P 

IV  THE  IMPULSE  RESPONSE  FUNCTIONS 

Using  the  results  of  section  III,  the  impulse  response  functions 

can  be  obtained  as  the  inverse  Fourier  transforms  of  the  frequency 

response  functions  u  (f)  and  u  '  (f) .  Decomposing  the  right  side  of  (16) 

P  P 

into  partial  fractions  gives 


u  (f)  = 


Ar 


Br 


cos2,n'f  -  z2 


cos2iTf  -  zT 


(26) 


where 


A  = 


(p/6)  Zp  +  (p/3) 

z+  -  Z~ 


and  Br 


(p/6)  Zp  +  (p/3) 

z~  —  z+ 
p  p 


(27) 


Then 


u  (t)  =  A  !  exp(2TTift)  df  +  B  j  exP  (2irif  t) 


COs2fTf  -  z" 


P  t  cos2irf  -  zx 

d  p 


df 


(28) 


P 

Tri 


z^dz 


z2-2z+  +  1 

P 


+  !p 

ttT 


z^dz 


oo 

E 


5 (t-i) 


^  z  -2z  z  +  1 

<*  p 

C  * 


(29) 


the  unit  circle,  C,  being  traversed  in  the  positive  sense  for  t>0.  To 
see  that  (29)  follows  from  (28)  consider  (12)  and  the  fact  that  the 
series  of  delta  functions  used  to  repeat  the  principal  part  of  u^(f) 
by  convolution  is  the  transform  of  the  series  in  t. 
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The  integrals  in  (29)  can  be  evaluated  by  the  method  of  residues  Cl] . 

The  roots  of  the  denominator  of  the  first  integrand  are  z++  and  z+~ , 

P  P 

and  for  the  second  z  and  z  ,  which  are  respectively  complex  conjugates 

P  P 

of  the  first  (reciprocal)  pair.  Then,  for  all  p>0,  the  impulse  response 
function 


u  (t) 

P 


(2K  z++|t|+  2 
P  P 


V 


-+  t 


CO 

)  .XI  <5  (t-i)  = 


l=-oo 


V  (t) 

P 


.£ 

l=_oc 


6 (t-i) , 


where 


K  = 


++ 

A  /  (z 
P  P 


++ 

Vzp  )  , 


and 


L„  = 


B  /(z 
P  P 


-+ 


Vzp  )  . 


(30) 


(31) 


(32) 


Using  the  same  conjugacy  relationships  between  the  roots,  for 
0<p<72 , 

Vp(t)  =  4Re{Kpz++ I t I }  =  4|Kp| |Zp+| I tlcos  (Arg  Kp  +  |t|Arg  Zp+)  (33) 

In  the  limit  as  p-K),  vp(t)  ->  0  uniformly  everywhere. 

For  p>72  there  are  four  distinct  negative  roots  occurring  in  reciprocal 
pairs.  The  impulse  response  function  becomes  the  sum  of  two  complex  numbers 
with  different  amplitudes  (decaying  exponentially  in  |t|)  rotating  in 
opposite  senses  around  the  t  axis,  piercing  the  imaginary-t  plane  at  half 
integral  values  of  t  and  the  real-t  plane  at  integral  values  of  t.  In  the 
limit  as  p-*00,  v  (t)  1  for  t=0  and  -*■  0  elsewhere.  That  is,  u  (t)  is  a 

ir  ir 

spike  at  the  origin. 

The  function  v^  (t)  must  be  evaluated  by  a  limit  process.  The  limits 
from  the  right  and  the  left  are  the  same  and  thus  Vp(t)  is  continuous  at 
p=72  for  all  values  of  t.  The  function  v  (0)  in  particular  is  plotted  in 
figure  4 . 
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The  functions  v^(t)  for  log  p  =  -3, -4, -5 

The  analysis  for  u  '  proceeds  as  for  u 

P  p 

and  B  '  are  different. 

P 


are  plotted  in  figure  5 . 
except  that  the  constants 


A  '  =  (2z+  -  2) / (z+  -  z") 
P  P  P  P 


and 


(34) 


-3 


V  =  (2zp '  2,/(zp "  ZP  (35) 

v  ' (0)  is  plotted  in  figure  4,  and  the  functions  v  ' (t)  for  log  p  = 
P  P 

-4,-5  are  plotted  in  figure  6. 


V  DISCUSSION  AND  AN  EXAMPLE 

Transforming  (15)  back  into  the  time  domain  and  reverting  to 
the  discrete  vector  symbolism  gives 

a  =  u*y-u'*  (y '  -  —  c'  ) ,  p  >  0,  (36) 

1  p  7  p  p  ^ 

and  the  operation  of  the  smoothing  spline  is  seen  to  be  representable 
in  the  time  domain  as  the  sum  of  two  convolutions.  The  first  term 
accurately  represents  the  operation  when  end  effects  can  be  ignored; 
for  instance  far  enough  away  from  the  ends  of  the  data  so  that  the 
filter  cannot  feel  them,  or  in  the  cases  in  which  a  reasonable  extra¬ 
polation  of  the  data  can  be  made.  The  second  term  is  peculiar  to  the 
smoothing  spline  and  accounts  for  its  apparent  lack  of  end  effects. 

It  is  necessary  to  consider  the  sense  in  which  the  smoothing  spline  has 
no  end  effects. 

This  phrase  is  appropriate  if  one  accepts  the  global  variational 
definition  of  the  smoothing  spline  as  implying  that  its  quasi-local 
behavior  is  the  same  at  the  ends  of  the  data  set  as  near  the  center. 
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Indeed  the  Euler-Lagrange  equations  must  be  satisfied  at  all  interior 


points.  In  this  sense  the  second  term  in  (36)  "fixes"  the  computations 

so  that  u  alone  is  an  adequate  description  of  its  operation,  in  which 
P 

case  one  need  only  consult  figure  1  in  selecting  an  appropriate  p  value. 
We  have  been  using  the  spline  in  this  sense.  It  is  more  appropriate  in 
this  paper,  however,  to  consider  the  end  behavior  in  terms  of  the  con¬ 
volution  representation. 

If  a  different  choice  had  been  made  in  extending  the  vectors  c,  y 
and  a,  the  contribution  of  the  second  term  in  (36)  would  have  been  dif¬ 
ferent.  In  fact  there  are  (an  infinite  number  of)  extensions  such  that 
this  second  term  would  contribute  nothing.  They  can  be  computed  by  con¬ 
volving  Up  with  y  and  requiring  that  the  resulting  values  be  equal  to  a^ 
for  0<i<n.  If  one  accepts  these  extensions  as  "reasonable",  then  again 
u^  by  itself  adequately  describes  the  filter.  This  implies  that  the 
smoothing  spline  can  be  used  as  a  predictor  or  to  achieve  higher  reso¬ 
lution  in  spectral  estimates.  If  the  physical  model  entails  that  the 
data  vary  according  to  some  principle  of  minimum  curvature  at  least  at 
lower  frequencies  such  an  idea  is  certainly  tenable,  but  as  in  other 
kinds  of  prediction  and  spectral  enhancement  (maximum  entropy,  for 
instance)  due  consideration  must  be  given  to  the  underlying  model. 
Prediction  is  reliable  only  if  the  model  is  appropriate. 

Four  smoothings  of  tree-ring  width  sequences  are  shown  in  figure  7. 
The  same  figure  appears  in  [2] .  The  dashed  lines  show  orthogonal  poly¬ 
nomial  fits  generated  by  the  INDXA  program  C53  the  orders  of  which  are 
determined  by  stopping  when  an  additional  order  fails  to  reduce  the 
residual  variance  by  5%  or  more.  The  log  p  value  was  -4  in  each  case. 
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Tree-ring  width  sequences  from  forest  interior  sites  tend  to  exhibit 
behavior  similar  to  that  seen  in  records  of  complex  physical  processes 
such  as  level  changes  in  streams  and  reservoirs  and  climatic  variations. 
These  data  exhibit  sudden  changes  in  level  (episodic  behavior)  over  a 
wide  range  of  time  scales.  The  spectra  tend  to  be  lumpy  and  have  concen¬ 
trations  of  power  near  zero  frequency.  The  particular  suitability  of  the 
smoothing  spline  as  a  low  pass  filter  for  certain  tree-ring  records,  and 
perhaps  for  other  similar  data,  is  accounted  for  by  the  shape  of  the 
impulse  response  functions.  The  spikey  exponential  envelope  with  the  |t| 
dependence  ensures  that  the  filter  will  be  responsive  to  sudden  changes 
in  level.  On  the  other  hand,  the  fact  that  it  is  a  low-pass  filter  with 
a  minimum  curvature  property  ensures  that  power  near  zero  frequency  will 
be  preserved.  The  frequency  response  function  is  smooth  at  the  origin. 
Thus  trends  and  slow  wandering  in  the  data  will  be  preserved.  Selecting 
a  value  of  S  (or  p)  affects  a  compromise  between  fitting  and  smoothing 
the  data.  We  see  how  this  compromise  manifests  itself  in  the  digital 
filters . 

Programming  of  equations  (31)  and  (4')  can  be  done  in  various  ways. 

We  have  not  tried  to  optimize  the  computations  beyond  using  an  equation 
solver  specialized  for  symmetric  band  matrices.  The  user  concerned  about 
storage  might  try  a  Toeplitz-specif ic  algorithm.  The  only  caution  regards 
the  possible  necessity  of  using  double  precision  accumulation  variables. 

We  found  that  in  single  precision  (32  binary  bits)  the  computations  were 

-4 

unstable  for  values  of  p  smaller  than  10  roughly. 
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APPENDIX  A 


It  is  not  obvious  that  the  weight  function  shown  in  figure  6  will, 
in  the  limit  of  small  p,  give  a  series  of  values  which  lie  on  a  straight 
line  for  O^i^n  and  be  zero  outside  this  interval  when  convolved  with 
anything  at  all.  The  purpose  of  this  appendix  is  to  examine  this  limit¬ 
ing  case. 

For  0<p<72,  and  thus  in  the  limit  as  p  -*■  0 , 

II  1 1 1 

v  ' (t)  =  4Re{Kp ' z^+ ' 1 }=4 | Kp ' | | z++ |  cos (Arg  Kp+|t|Arg  z^+) .  (Al) 

Letting  r,  0'  and  0  denote  small  positive  numbers,  for  small  p  Al  can 
be  written 

v  ' (t)  =  4  Ik  '  |  (l-r)  1 1 1  cos  (  (3tt/4)  -0'  -  1 1 1  0 )  .  (A2) 

p  P 

Expanding  in  r  and  0 '  and  0 , 

v  1  (t)  =  k  [1  -  1 1 1  r  +  o(r2)]  El-0'  -  1 1 1  0  -  o(02)],  (A3) 

P 

where 

k  =  -4|K  '  |/V2. 

To  first  order  therefore 

v  '  (t)  =  k(l-0‘-  s  1 1 1 )  ,  for  s  <  <  1 1 1  <  <  1/s,  (A4) 

ir 


where  s  =  r  +  0. 


1 


Denoting  the  four  non-zero  values  of  (y'-vc  ')  by  w  , 

ir  P  “  1 

wn+l'  at  least  the  following  four  requirements  must  be  met: 

k(l-0,)w_1  +  k (1-0 ' -s) wQ  =  0, 

k(l-@'-s)w_1  +  k (1-0 ' ) wQ  =  aQ, 

k(l-0')wn  +  k(l-0'-s)wn+1  =  an, 

and  k(l-0'-s)w  +  k(l-0')w  =  0. 

n  n+1 


w 


O' 


w 


n 


(A5) 

(A6) 

(A7) 

(A8) 
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Since  in  the  limit  of  small  p,|K  1 |  and  thus  |k|  become  very  large, 

ir 

k  can  be  divided  out  of  A5  -  A8.  Solving  the  resulting  equations  for 

w  and  w  (A5  and  A6)  and  w  and  w  (A7  and  A8)  yields  to  first  order 
-1  0  n  n+1 

in  s , 


w_x  =  -aQ (1-0 ' -s) /2s ( 1-0 ' ) , 

(A9) 

wQ  =  a0/2s. 

(A10) 

w  -  a  /2s, 
n  nx  ' 

(All) 

and 

wn+l  =  -an(l-0,-s)/2s(l-0') . 

(A12) 

Since  s  >  0  for  all  finite  p  these  solutions  are  defined  for  the  limit. 
Also,  by  the  symmetric  linear  (to  first  order)  dependence  of  v  ' (t) 

IT 

on  t,  it  follows  from  A5  that  a^  =  0  for  all  -1/s  <<i<  -1;  from  A8, 
that  a^  =  0  for  all  n+l<i<<l/s;  and  from  A6  and  A7  that  a^  varies 
linearly  with  i  (to  first  order)  for  O^i^n.  In  the  limit  as  p-*0  the 
first  order  expressions  become  exact  and  l/s-*-°°. 
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APPENDIX  B 


The  first  terms  in  the  central  finite  difference  expansions  of 
the  second  and  fourth  derivative  operators  are  obvious  in  the  matrices 
Q  and  R.  Corresponding  to  equations  (3')  and  (4')  there  are  the  coupled 
differential  equations 

CD4  +  (p/3)  (D2  +  6)]  C(t)  =  pD2Y(t)  (Bl) 

and 

A  ( t)  =  y  -  —  D2C  ( t)  ,  p>0 .  (B2) 

P 

Conversely,  a  finite  difference  solution  of  (Bl)  and  (B2)  will 
satisfy  equations  (3')  and  (4').  More  generally,  the  theory  of  Reinsch's 
smoothing  splines  is  related  to  the  theory  of  the  interpolated  finite 
difference  solutions  of  certain  ordinary  differential  equations  with 
constant  coefficients.  The  remainder  of  this  appendix  is  a  brief 
analysis  of  the  equations  (Bl)  and  (B2) ,  analogous  to  that  done  for 
the  spline  equations,  which  will  facilitate  some  comparisons  the  reader 
may  wish  to  make  between  the  continuous  and  discrete  systems. 

Fourier  transforming  equations  (Bl)  and  (B2)  and  solving  for  A(f) 
in  terms  of  Y(f) ,  gives 


p(3-27T2f2) 

A  (f )  =  - — - — —  Y(f)  =U  (f)Y(f).  (B3) 

24irf  +p(3-2tt  f2) 

Note  that  this  same  transform  would  arise  from  the  single  equation 


(D4  +  (p/3)  (D2  +  6))A(t)  =  (p/3)  (P2  +  6)  Y  ( t)  . 


(B4) 


This  won't  work  in  the  discrete  case  because  of  the  correction  terms. 
The  roots  of  the  denominator  for  0$p-$72  are 


fp+  =  (P  +  i  (p  (72-p)  )  1/2 )  /24tt2 

cin  d  i  ^  p 

f2_  =  (p  -  i  (p  (72-p)  )  )/24tt 

ir 


(B5) 

(B6) 
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and  therefore 


and 


+<£p+)1/2- 


-if2*) 1/2 

P 


-(f2_)1/2. 

P 


(B7) 

(B8) 

(B9) 


(BIO) 


Corresponding  to  the  constants  A  and  B 

P  P 


there  are 


2  2+  4  9+  o_ 

a  =  p  (3-2tt  f  )  /24tt  ( r  +  -  f  ), 
P  P  P  P 


and 


2  2-  4  2-  2+ 

3  =  p(3-2rr  f  )  /24ir  (f  -  f  )  , 

P  P  ir  P 


in  the  inverse  transform 


(Bll) 


(B12) 


while 


k 

P 


Vf 


++ 
P  ' 


(B13) 


and 


£P  -  vf;+ 


(B14) 


The  contour  of  integration  consists  of  the  real  axis  and  the  upper 
half  of  the  point  at  infinity,  traversed  counterclockwise  for  t>0. 

For  0<p<72 

U  (t)  =  -2iTexp  (-2iTlm{  f^+}  1 1 1 )  Im{k^exp  (27riRe{fp+}  1 1 1 )  }  .  (B15) 

For  p^72,  each  contour  (t<0)  includes  two  of  the  real  roots  and 

U  (t)  =  -ir(k  sin  2tt  f 1 1 1  +  1  sin  2ir  f_+ 1 1 1 )  .  (B16) 

P  P  P  P  P 

There  is  no  decay  in  the  continuous  case  for  p^72. 

Evaluation  of  the  limits  as  p  approaches  72  from  the  left  and  right 

will  show  that  for  any  value  of  t,  U  (t)  is  continuous  at  p=72. 

P 


18 


ACKNOWLEDGEMENTS 


The  authors  wish  to  thank  Paul  Stoffa,  John  Shepherd  and  Keith 
McCamy  for  many  constructive  criticisms.  Funds  for  this  work  were  pro¬ 
vided  through  the  Climate  Dynamics  Program  of  the  National  Science  Founda¬ 
tion,  grant  no.  ATM-77-19217. 


19 


FIGURE  CAPTIONS 


Figure  1. 


Figure  2. 


Figure  3. 


Figure  4. 


Figure  5. 


Figure  6 . 


Figure  7. 


The  gain  of  the  frequency  response  function  u  (f)  for  various 

P 

values  of  the  parameter  p  over  the  range  .001^f^.5. 

The  base  ten  logarithm  of  the  negative  of  the  frequency  response 
function  u  ' (f)  for  various  values  of  the  parameter  p  over  the 

.c' 

range  . 001^f<.5.  The  dashed  line  passes  through  the  maximi 
(circles)  on  these  curves. 

The  root  z++  of  the  characteristic  polynomial  r  (z)  .  Values  of 
P  P 

the  parameter  p  are  shown  along  the  curve . 

The  central  values  vp(0)  and  v^ 1  (0)  of  the  smooth  parts  of  the 
impulse  response  functions  u^  and  u ' . 

The  smooth  part  v^(t)  of  the  impulse  response  function  u^  for 
three  values  of  p  most  useful  for  standardizing  forest  interior 
tree-ring  series. 

The  smooth  part  v  ' (t)  of  the  impulse  response  function  u^ ' 
for  three  values  of  p  most  useful  for  standardizing  forest 
interior  tree-ring  series. 

Four  examples  of  the  smoothing  spline  applied  to  forest  interior 
ring-width  series.  Each  spline  was  computed  for  log  p  =  -4.0. 
The  solud  curve  is  the  spline  fit  and  the  dashed  curve  is  the 
orthogonal  polynomial  fit  from  program  INDXA  (see  ref.  4) .  The 
four  series  are  from  the  same  site  in  southeastern  New  York  and 
the  lower  two  plots  are  curves  from  the  same  tree. 
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